!
!     Kombine     Combine the set of frames.
!     Copyright (C) 2006  Filip Hroch, Masaryk University, Brno, CZ
!
! $Id$
!

Program SKombine

!
! Experimental version, subpixel transformation
!
! Sum (combine) the frames.
! To use generate .TRA file as the output of the munimatch (with -t switch)
!
!

  use interpol
  use optfile
  use stat
  use daoio

  implicit none

! Compiled-in parameters (perhaps, no need to adjust)

! Version:
  Character(len=*), parameter :: version = 'SKOMBINE Version 0.4, (C)1997-08 F.Hroch, Masaryk University,Brno,CZ '

! Minimum of values
  Real :: minvalue = 0.0

! filename and line lenght parameters
  integer, parameter :: lenfilename = 1024, lenline = 80  

! max stars to intensity computing
  integer, parameter :: maxstar = 256

! FITS
  Logical :: simple,extend,anyf,exist
  Integer :: n, nalldate, maxalldate
  Integer :: i, j, ii, istat, bpix, maxvalue, bitpix = -32, naxis, naxes(2) 
  Integer :: BlockSize,pcount,gcount
  character(len=8) :: FILTER='FILTER',DATE_OBS='DATE-OBS',EXPTIME='EXPTIME', &
       OBJECT='OBJECT'

  Real, dimension(:,:), allocatable :: KOMBINE, IMAGE
  real, dimension(:), allocatable :: obr
  Real :: bscale, bzero
  Real :: oneexp, etime
  Real :: xx,yy,c,r(2),m(2,2),t1(2),t2(2),x,y,f,fs,sky,amag
  real, dimension(maxstar) :: fx, fx0, asky, asky0
  integer :: ix, iy, nl
  character(len=lenline) :: onedate, obj, filtr
  character(len=lenline), dimension(:), allocatable :: alldate, alld
  character(len=lenfilename) :: sb1,sb2,sb,sbtra,sbmat
  character(len=lenline) :: buf,key,value
  logical :: initialize = .true.

! selection of method of interpolation
! 0=no interpolation, 1=bilinear, 2=bicubic, 3=bspline
!   fast                slow        slower     hung up
!   tooths              good        perfect    little bit over perfect
  integer :: sint = 0

! set rescale factor (rescale*dimensions of first image, if nx,ny < 0)
  real :: rescale = 1

! set dimensions of output image 
  Integer ::  nx = -1, ny = -1
  
! Output image name:
  Character(len=lenfilename) :: namekombi='skombine.fits'

! Flux recompute to one level
  logical :: flux = .false.

! Normalize output,
! 0.. no normalization, 1 .. mean (sum divided by number), 2.. normlize to sky
  integer :: norm = 0

! Time start of first exposure
  character(len=lenline) :: first = ''

!------

!  write(*,*) trim(version)

! init variables
  istat = 0
  obj = ''
  etime = 0.0
  nalldate = 0
  maxalldate = 2
  f = 1.0
  fs = 0.0
  sky = 0.0
  allocate(alldate(maxalldate))

! check options file
  if( opt_open(21,'kombine.opt',i) )then
     do
        call opt_read(21,key,value,i)
        if( i == opt_end .or. i == opt_err ) exit
        !write(*,*) trim(key),':',trim(value)

        if( key == 'Method of interpolation' ) then
           read(value,*,iostat=i) sint
        endif

        if( key == 'Rescale factor' ) then
           read(value,*,iostat=i) rescale
        endif
        
        if( key == 'Output normalized by' ) then
           read(value,*,iostat=i) norm
        endif

        if( key == 'Multiply flux' ) then
           read(value,*,iostat=i) flux
        endif
        
        if( key == 'Output width' ) then
           read(value,*,iostat=i) nx
        endif

        if( key == 'Output height' ) then
           read(value,*,iostat=i) ny
        endif

        if( key == 'Output BITPIX' ) then
           read(value,*,iostat=i) bitpix
        endif

        if( key == 'Output filename' ) then
           namekombi = value
        endif

        if( key == 'Key DATE-OBS' ) then
           read(value,*,iostat=i) DATE_OBS
        endif

        if( key == 'Key EXPTIME' ) then
           read(value,*,iostat=i) EXPTIME
        endif

        if( key == 'Key FILTER' ) then
           read(value,*,iostat=i) FILTER
        endif

        if( key == 'Key OBJECT' ) then
           read(value,*,iostat=i) OBJECT
        endif

     enddo
     call opt_close(21)
  end if

  ! input of top line with init values
  read(*,*,end=666,iostat=j) sb,x,bpix,n,i,ii
  if( j /= 0 ) stop 'Bad format of input data.'
  if( sb /= '' ) namekombi = sb
  if( x > epsilon(0.0) ) rescale = x
  if( bpix /= 0 ) bitpix = bpix
  if( n >= 0 ) norm = n
  if( i >= 0 ) sint = i
  if( ii > 0 ) flux = .true.
  
  n = 0
  istat = 0
! input filename
  Do
     read(*,*,end=20,iostat=i) sb
     if( i /= 0 ) stop 'Bad format of input filename.'
     write(*,*) trim(sb),':'
         
! input initialisation of the current file
     Call ftopen(25,sb,0,BlockSize,istat)
     Call ftghpr(25,2,simple,bpix,naxis,naxes,pcount,gcount,extend,istat)
     if( istat /= 0 ) goto 666
     CALL ftgkye(25,'BSCALE',bscale,buf,istat)
     CALL ftgkye(25,'BZERO',bzero,buf,istat)
     !        BSCALE, BZERO not presented
     if( istat /= 0 ) then
        istat = 0
        bscale = 1.0
        bzero = 0.0
     endif
     call ftgkye(25,EXPTIME,oneexp, buf, istat)
     if( istat /= 0 ) then
        oneexp = -1.0
        istat = 0
     endif
     call ftgkys(25,DATE_OBS,onedate, buf, istat)
     if( istat /= 0 ) then
        onedate = ''
        istat = 0
     endif
     call ftgkys(25,OBJECT,obj, buf, istat)
     if( istat /= 0 ) then
        obj = ''
        istat = 0
     endif
     call ftgkys(25,FILTER,filtr, buf, istat)
     if( istat /= 0 ) then
        filtr = ''
        istat = 0
     endif

! image - data input
     allocate(IMAGE(naxes(1),naxes(2)))
     !IMAGE = 0.0
     call FTG2DE(25,1,minvalue,naxes(1),naxes(1),naxes(2),IMAGE,anyf,istat)
     if( istat /= 0 ) goto 666
     call ftclos(25,istat)
     if( istat /= 0 ) goto 666

! image - transformation input
     sbtra = sb
     call another(sbtra,'.fits','.TRA')
     open(27,file=sbtra,status='old',iostat=i)
     if( i == 0 ) then
        read(27,*,iostat=i) sb1,sb2,c,r,m,t1,t2
        close(27)
     endif
     if( i /= 0 ) then
        write(*,*) 'No such TRA file or invalid format. Skipped.'
        goto 665
     endif

! image - determine flux parameters
     if( flux ) then
        
        sbmat = sb
        call another(sbmat,'.fits','.MAT')
        j = 0
        call dao_rdhead(28,sbmat,nl,j)
        i = 0
        do
           call dao_rdstar(28, nl, ii, x, y, amag, xx, j)
           if( j == dao_end ) exit
           !write(*,*) i, x, y, amag(i), asky(i)
           if( j == dao_ok ) then
              i = i + 1
              fx(i) = 10**(-0.4*(amag - 25.0))
              asky(i) = xx
           endif
           if( i == maxstar ) exit
        enddo
        call dao_close(28,j)
        if( i > 0 .and. .not. initialize ) then
           where( fx > 0 )
              fx = fx0/fx
           end where
           
           call rmean(fx(1:i),f,fs)
           asky = f*asky - asky0
           call rmean(asky(1:i),sky,x)
           !write(*,*) f,sky
        else
           f = 1.0
           fs = 0.0
           sky = 0.0
        end if
     end if

! first image - internal buffers initialisation
     if( initialize ) then
        initialize = .false.

        if( nx < 0 .and. ny < 0 )then
           nx = nint(rescale*naxes(1))
           ny = nint(rescale*naxes(2))
        endif

        allocate(KOMBINE(nx,ny))
        KOMBINE = 0.0

        if( flux ) then
           fx0 = fx
           asky0 = asky
           f = 1.0
           fs = 0.0
           sky = 0.0
        end if

     end if

     if( sint == 3 ) then
        ! init bspline
        call bspline_init(IMAGE)
     endif

     ! compute sumation
     do i = 1, nx
        do j = 1, ny
           
           xx = i/rescale - t1(1)
           yy = j/rescale - t1(2)
           x = c*(r(1) + xx*m(1,1) + yy*m(1,2)) + t2(1)
           y = c*(r(2) + xx*m(2,1) + yy*m(2,2)) + t2(2)

           !write(*,*) i,j,x,y
              
           select case (sint)

           case(0)
              ! no interpolation
              ix = nint(x)
              iy = nint(y)
              if(  1 <= ix .and. ix <= size(IMAGE,1) .and. &
                   1 <= iy .and. iy <= size(IMAGE,2) )then
                 KOMBINE(i,j) = KOMBINE(i,j) + f*IMAGE(ix,iy)-sky
              end if
                         
           case(1)
              ! bilinear interpolation
              KOMBINE(i,j) = KOMBINE(i,j) + f*bilinear(x,y,IMAGE)-sky
              
           case(2)
              ! bicubic interpolation
              KOMBINE(i,j) = KOMBINE(i,j) + f*bicubic(x,y,IMAGE)-sky
              
           case(3)
              ! bspline interpolation
              KOMBINE(i,j) = KOMBINE(i,j) + f*bspline_func(x,y,IMAGE)-sky
              
           end select

        enddo
     enddo
     if( sint == 3 ) then
        call bspline_shutdown
     endif

     if( onedate /= '' .and. oneexp > 0.0 ) then
        write(buf,'(a,1x,f0.2,1x,a,1x,f0.3,a,f0.3,1x,f0.2)',iostat=i) &
                   trim(onedate),oneexp,trim(filtr),f,'+-',fs,sky
        if( i /= 0 ) buf = trim(onedate)//' ..'
        etime = etime + oneexp
        if( first == '' ) first = onedate
        
        nalldate = nalldate + 1
        if( nalldate > maxalldate ) then
           allocate(alld(maxalldate))
           alld = alldate
           deallocate(alldate)
           maxalldate = 2*maxalldate
           allocate(alldate(maxalldate))
           alldate(1:maxalldate/2) = alld
           deallocate(alld)
        endif
        alldate(nalldate) = buf 
        
     endif
     n = n + 1
665  continue
     deallocate(IMAGE)
  enddo
20 continue

  if( n == 0 ) stop 'No input image.'
  if( bitpix == 0 ) bitpix = bpix

  ! normalize
  select case (norm)
  case(1)
     KOMBINE = KOMBINE/n
  case(2)
     allocate(obr(nx*ny))
     obr = pack(KOMBINE,.true.)
     call rmean(obr,sky,y)
     KOMBINE = KOMBINE / sky
     deallocate(obr)
  end select

  ! range checking
  if( bitpix > 0 )then
     maxvalue = 2**bitpix - 1
     where( KOMBINE < minvalue )
        KOMBINE = minvalue
     endwhere
     where( KOMBINE > maxvalue )
        KOMBINE = maxvalue
     endwhere
  end if

! form the output image
  Inquire(file=namekombi,exist=exist)
  if( exist ) then
     open(26,file=namekombi,iostat=i)
     close(26,status='DELETE', iostat=i)
     if( i /= 0 ) then
        write(*,*)"Error occured during remove of output file."
        goto 666
     endif
  endif
  Call ftinit(26,namekombi,1,istat)
  naxis = 2
  naxes = (/nx, ny/)
  Call ftphps(26,bitpix,naxis,naxes,istat)
  if( bitpix > 0 ) then
     call ftpkye(26,'BSCALE',1.0,10,'',istat)
     call ftpkye(26,'BZERO',2.0**(bitpix-1),10,'',istat)
  endif
  Call ftpkys(26,'OBJECT',obj,'',istat)
  Call ftpkys(26,'FILTER',filtr,'',istat)
  call ftpkye(26,'EXPTIME',etime,10,'Sum of all times',istat)
  if( first /= '' ) call ftpkys(26,'DATE-OBS',first,'UTC of serie exposure start',istat)
  if( nalldate > 0 ) then     
     write(buf,'(a,i0,a)',iostat=i) &
          'The image is a result of composition of ',nalldate,' exposure(s).'
     if( i== 0 ) call ftpcom(26,buf,istat) 
     call ftpcom(26,'The time of start, the exposure time, the filter,',istat) 
     call ftpcom(26,'the intensity multiplier and relative sky value for each image used:',istat)
     do i = 1, nalldate
        call ftpcom(26,alldate(i),istat)
     enddo
     call ftpcom(26,'All times are derived from DATE-OBS and EXPTIME keywords.',istat)
  endif
  if( norm > 0 ) then
     select case (norm)
     case(1)
        buf = 'Fluxes are normalized by number of images.'
        i = 0
     case(2)
        write(buf,'(a,f0.3)',iostat=i)'Fluxes are normalized by sky value: ',sky
     end select
     if( i == 0 ) call ftpcom(26,buf,istat)
  end if
  select case (sint)
  case(0)
     buf = 'No interpolation has been used.'
  case(1)
     buf = 'Bilinear interpolation has been used.'
  case(2)
     buf = 'Bicubic interpolation has been used.'
  case(3)
     buf = 'Bspline interpolation has been used.'
  end select
  call ftpcom(26,buf,istat)
  call ftpcom(26,version,istat)
  if( istat /= 0 ) goto 666
  
  Call ftp2de(26,1,nx,nx,ny,KOMBINE,istat)
  
  CALL ftclos(26,istat)
  if( istat /= 0 ) goto 666

  ! print output info
  write(*,*) 'Output image: ',trim(namekombi)
  write(*,*) 'Dimension: ',nx,'x',ny

666 continue
  if( allocated(KOMBINE) ) deallocate(KOMBINE)

  if( istat /= 0 )then
     write(*,*) 'The following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif

End Program SKombine

